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Foreword 


A need to investigate solutions to the inhomogeneous equation of heat conduc- 
tion with time-dependent sources and the so-called linearized radiation (or 
insulation) boundary condition arose directly from a study of the constraints on 
lunar thermal history posed by systematic analysis of returned lunar samples and 
geophysical data reported by Conel et al. in 1972 (see Ref. 1). A great deal of 
numerical modelling was carried out to support and guide us to the principal 
conclusions presented in that report; these specific model results will be published 
elsewhere. The present document describes only the mathematical problem 
involved and its numerical solution. 

Studies of planetary thermal history have, with time, evolved mathematical 
models of ever-increasing complexity. Added complications have involved inclu- 
sion of radiative transfer with the (radiative) thermal conductivity as a prescribed 
function of temperature, inclusion of latent heat associated with phase changes, 
simulated convection in molten and solid zones, redistribution of heat sources 
upon solidification according to specified laws, etc. Most of these refinements 
render the mathematical problem nonlinear, and hence force numerical integra- 
tion of the heat equation from the outset. Sophisticated models also require 
specification of an increasing number of material parameters that are often poorly 
determined. In my opinion, it is hard to justify an increase in model complexity 
when the fundamental data to be used are not well known. An increase in the 
degrees of freedom accompanying construction of elaborate model schemes also 
allows an investigator to achieve his desired result by more and more diverse 
means. To avoid some of these possible sources of difficulty, I have purposely 
chosen to deal with comparatively simple analytic solutions to the equation of 
heat conduction. My attitude is that until it is compellingly shown that elementary 
procedures and simple assumptions fail to explain the observations, it is not 
worthwhile abandoning such procedures and assumptions. 

Ordinarily, there is no place for publication of the details of complex computer 
codes in the scientific literature. I consider this an unfortunate necessity, since it 
may mean that there is no way of checking the procedure or accuracy of a lengthy 
numerical exercise, and thus no basis for judging its validity. I have prepared the 
present report to compensate in part for any such deficiency in my own work. X 
hope in addition that the documentation and code may prove useful to others 
with their special problems. 
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Abstract 


A computer program (Program SPHERE) solving the inhomogeneous equation 
of heat conduction with radiation boundary condition on a thermally homoge- 
neous sphere is described* The source terms are taken to be exponential functions 
of the time. Thermal properties are independent of temperature. The solutions 
are appropriate to studying certain classes of planetary thermal history. Special 
application to the Moon is discussed. 
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Solution of the Equation of Heat 
Conduction with Time-Dependent 
Sources: Programmed Application 
to Planetary Thermal History 


L Introduction 

This report is a concise documentation of a programmed solution to the inhomo- 
geneous equation of heat conduction in spherical coordinates with time-dependent 
sources (Program SPHERE). The application here is in the study of planetary 
thermal history, and in particular the Moon (see Ref. 1), The sphere is homo- 
geneous in density and thermal properties . The particular example given is 
specialized from a general case for radial distribution of heat sources and initial 
temperature given by Lowan (Ref. 2). The boundary condition at the outer surface 
is the so-called linearized radiation boundary condition or “Fourier's problem of 
the third kind”; in this particular instance the body is considered to radiate to a 
medium at constant temperature. This condition also applies to problems where 
a thin skin of insulating material exists on the exterior of the body (Ref. 3) such 
that the heat capacity of the skin can be neglected. This is tantamount to saying 
that if a change in temperature occurs on the inner insulating boundary, then 
the exterior medium responds "instantaneously” to establish a linear temperature 
gradient in the insulation itself. In this case the outer surface temperature of the 
insulation is held fixed, which corresponds mathematically to the case previously 
discussed where radiation is to a medium at constant temperature. 

Two points should be emphasized in applying these results: (1) The surface 
temperature solved for is the temperature of the medium just beneath the 
insulating layer, and (2) Instead of dealing with the more complicated problem 
of varying surface temperature (i.e., the temperature at the outer physical surface), 
we have chosen the boundary temperature to be fixed, in most cases taking it to 
be the average value of whatever the expected sinusoidal variation might be on 
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an atmosphereless body rotating with given angular velocity about its own axis, 
and around the Sun. The boundary temperature, of course, can be assigned 
arbitrarily. Thus, the details of the fourth power nonlinearity in the usual bound- 
ary conditions are not dealt with; at the same time, these details are completely 
unimportant in understanding thermal problems of deep planetary interiors. (They 
may, of course, become important in applications where near-surface regolith 
conditions are of interest; in this case numerical techniques would be advisable 
from the outset, although viable alternative procedures, still requiring computers 
for hard answers, are useful (see Ref. 4 for details)). 

The reason for utilizing this general approach in analyzing problems involving 
insulation is that the full, more complicated problem of dealing with the thermal 
properties of the outer layer is avoided by a simple mathematical trick. At the 
same time, we have been forced to neglect any radioactively generated heat 
contribution in the insulating blanket itself. This could normally be dealt with in 
practice by an approximate separate calculation of equating sources to flow, since 
the insulating layer is usually thin relative to the planet ary radius. 

An additional feature has been added to the original Lowan calculations. In 
considering the physical problem of lunar thermal history, it became evident that 
the boundary condition might necessarily be a function of time. The specific 
problem treated is that of having the insulation vary in thickness as a step function 
of time at some specified time f > 0. In this instance, t = 0 corresponds to the 
time of formation of the Moon, 4.6 X 10 9 years ago. In the lunar example, 
t>~ 1 X 10 9 years (3.1536 X 10 16 s), although in general the value of if may be 
anything, in many practical circumstances zero. When a nonzero if is used in the 
program, the new origin of time coordinates is such that t — 0 corresponds to t', 
and the calculations are made in specified, time incremental steps At with the 
first being f + At. Thus, for example, if t' = 1 X 10 9 years, the time corresponding 
to the present would be 3.6 X 10 9 years, and if 10 time-increments were needed. 
At — 3.6 X 10 8 years. The first calculated value would correspond to an actual 
time from the standpoint of radioactive source strength of 1.0 X 10 9 years plus 
3.6 X 10 s years or 1.36 X 10 9 years. Once t f and At have been assigned, a specific 
portion of the tiine-history may be examined by redefining the origin of time, 
and specifying N y the number of time steps, in a suitable fashion. Suppose, for 
example, that the thermal regime at a single time 4.5 X 10 9 years is wanted. The 
time origin (program statement 999) is simply redefined as TIME = 1.41912D17 — 
At, and N is set equal to one. 


II. Statement of Problem 

The formal boundary value problem for temperature T as a function of r and f 
which we solve, is as follows: 


dT K 2 3T K d z T _ 

r dr pc dr'- “ 


( 1 ) 


lim T(r,t) = f(r) (2) 
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( 3 ) 



*L + hT = hT s , f = R (0 <t<f) (4) 

In Eqs. (1) through (4), 

p = density in g/cm -3 
c = specific heat in cal g^ 1 °K“ 1 
K = thermal conductivity in cal cm -1 sec - ” 1 °K~ 1 

The factor h = K'/Kd , where 

K' ~ thermal conductivity of insulation 
d = thickness of insulation in cm 

At time f the thickness of insulation is assumed to change to d' y and Eq. (4) 
becomes 


+ h'T = h'T s , t > f (5) 

of 

The source function pertains to heat production from exponential decay of 
radioactive sources and is given by 

+(rj) = — y pA,(f ,r ) fty exp [ A|(f - *)] (6) 

where A;(T,r) are abundances (in g/g) of radionuclide species j at time T = 4,6 X 
10° years, i.e., the present. The Jty are rates of heat generation in cal g _1 sec -1 
and \j are the decay constants. The radionuclides considered are 40 K, 2S2 Th, 235 U, 
and 238 U (/ — 1, 2,3,4). The following relations between these have been assumed 
or are to be specified: 


F 232 Thl p 35 Ul 1 [Kl 

[«K] = 1.19 X 10- [K]; 1 T y F L = 4; ^ ; j jjj = P (7) 

The parameter may be specified arbitrarily but has not in the present problem 
been taken as a function of r. The final variable of the source function is [U], and 
may be specified arbitrarily. In all of the above, [■] signifies concentration in g/g. 


III. Solution of the Differential Equation 

As a result of the boundary condition (Eq. 5), the solution for T(r,t) subsequent 
to V requires the solution of a new boundary-value problem with a new "initial” 
temperature distribution f (r) appropriate to whatever problem is embodied in the 


JPL TECHNICAL MEMORANDUM 33-718 


3 



solutions to Eqs. (1) through (4) and Eq. (5). The formal solution to the system of 
Eqs. (1) through (4) is given by Lowan (Ref. 2; see his Eq. 2T) as 


00 — 


(h 2 + A») sin A „r 


T(r,t)-T,= — V _ 

r ^ R\i + h(Rh + 1) 

n = i 

+ j $ sin exp(-K\lt) dr| 


exp ( - k AS#) | if(%) sin A 


( 8 ) 


where 


h 


1_ 

Kd R 


and the kn are roots of 


(9) 


a COt aR + h — 0 (10) 

The conductivity of the insulating layer is K'; R is the planetary radius. It is 
ordinarily not known how many roots of this equation may be required to achieve 
a given stability to the convergence of the series in Eq. (8), so it is possible to 
specify the maximum number of solutions to Eq. (10) independently. These are 
ordinarily computed first and stored. In a typical problem using double-precision 
arithmetic, 75 to 100 terms may be summed for any given radial and time incre- 
ment, so that at least this many \ n are required. As a safety factor, the value 
MAXIT, defined in the Appendix, is usually given a value like 500. The summation 
over n is continued until ten consecutive values in the series give values differing 
by less than I0“ 3 . This criterion naturally depends upon uniform convergence of 
the series. The question of nonuniformity of convergence has not been formally 
investigated. Tkusfar, however, no numerical instability has been noted in any 
calculation carried out to date. 


IV. Specification of Initial Temperature and Distribution of 
Radioactivity 

In problems dealing with the Moon, we have considered spherical shell geome- 
tries like that shown in Fig. 1. 

For initial temperature, the Moon is divided into two radial, spherically sym- 
metric zones: 


f(i) = T(r,0) = T 0 + [T m ( ri ) - T c ] £ 0 < r < fi (11) 

= T m (R) + K(R - r)„ r x <r<R (12) 

Here, T 0 is equal to the initial temperature of accreting material (on the inde- 
pendent spherical accretion model) adjusted appropriately for any increase due 
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Fig. 1. Three-shell radioactive sphere 



Fig. 2. Initial temperature and radionuclide distributions 


to adiabatic compression. At 1 AU at the present time, for example, uncompressed 
material would have a temperature of about 400°K while, for lunar mass, adiabatic 
compression would account for about a 50 °K increase. The central temperature 
T 0 for such problems is thus near 450°K. In Eqs. (11) and (12), T m (r a ) is the 
melting point of postulated lunar material at radius r u T m (R) is the melting point 
at the surface, and T' m is the melting point gradient in °K/cm. 

The graphical form of Eqs. (11) and (12) describing independent accretion and 
surficial melting and differentiation is shown in Fig. 2. For numerical purposes, 
these quantities, where required, have been taken from experimental values 
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obtained from returned lunar samples. Thus, Ringwood and Essene (Ref. 5) would 
give the following values for Apollo 11 basalt: 


T m (R) = 1360° K 
TJr x ) - 1480°K 

T' m = 4.8 X 10- 6 °K/cm 


Equations (11) and (12) may be modified to describe other models of planetary 
origin. If a lunar-sized object is taken to be at the melting point throughout 
initially, the starting temperature in the absence of phase changes will be quite 
nearly a parabolic function of radius. To describe such a temperature profile, in 
Eq. (11) we set r x equal to R, and T m {r x ) equal to T m (R) y the melting point at the 
surface; in Eq. (12), T m (R) = 7^ = 0. The quantity T 0 is, then interpreted as the 
melting point at r = 0. 

The abundances of radioactive species on the accretion model are taken to be 
constant values within each zone, not exponentially decreasing functions of depth 
as assumed, for example, by Hanks and Anderson (Ref. 6). While there is field 
evidence for decrease of radioactive sources with depth in the Earth according 
to a law such as [U] — [U 0 ] (z positive downward), this has been substan- 
tiated only for crustal regions and most firmly in relatively localized batholithic 
portions of the crust. At any rate, if the portions of the Moon involved in any 
primordial melting are small, i.e., R — r 2 < < R ? then the detailed distribution of 
sources near the surface is of no consequence as far as the deeper temperature is 
concerned. 

The value of [Ui] is equal to the present "primordial” t7 abundance, [U 2 ] is 
equal to the present C7 abundance of the depleted zone, and [U 3 ] is equal to the 
present U abundance in the “crust,” and determined as follows: [UJ is the assumed 
primordial value of "untouched” or undifferentiated lunar material; [U 2 ] and 
[U 3 ] then follow by mass conservation arguments once it has been hypothesized 
what fraction / of [Ui] has been removed from the bleached zone by a differen- 
tiation process. If complete bleaching has occurred, then / — 1 for zone two and 
the entire mass of U in the zone between r 2 and t 1 has been concentrated into the 
zone between r 2 and R, 

With the concentration [UJ specified, together with the factor /, we get for a 
homogeneous body 

[U 2 ] = (1-/)[U 1 ] (13) 

and 


[UJ = 


R 3 “ f 2 + f( T l 

R*-rl 


Note that Eqs. (13) and (14) allow treatment of a variety of two-zone models as 
well. If the uranium abundance of postulated primordial material is [U a ] and the 
planet is conceived to differentiate into two portions with core of radius r 2 and 
(present) abundance [Ui], we set / equal to 1 — [Ui]/[UJ, and fi equal to 
zero. The factor / must be calculated by hand. It does not matter mathematically 
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that [UJ is specified for a region of zero volume. To obtain the temperature 
increase from radioactive sources alone, T s must be set equal to zero, and the 
initial temperature distribution taken zero throughout as well. 

Specification of the U abundance has been emphasized here and we have relied 
upon connections between U, Th, and K cited earlier to specify abundances of 
the other species involved, and their corresponding heat generations. It may well 
be that such systematic connections do not exist in all instances (except the 
235 U/ 238 U ratio). Thus /?, the [K]/[U] ratio from the orbital gamma ray data, 
can be shown to vary laterally over the Moon s surface. There is also no guarantee 
that it is constant throughout the lunar interior. A similar situation exists on 
the Earth where xenoliths from the deep interior (Ref. 7) have lunar-like 
rather than terrestrial-like [K]/[U] ratios. So the real distributions may be ones of 
great complexity. 

We have noted that there is room for considerable computational flexibility in 
the distributions of initial temperature or radioactivity in spite of specific analytic 
forms taken here. The user is cautioned, however, that calculations for decay of 
initial temperature or temperature changes for radioactive heating must be done 
sequentially unless boundaries of zones in the initial profiles coincide, as they 
are shown to do in Figs. 1 and 2. The assumption of spherical symmetry is always 
made. 

Integration of Eq. (8) is routine for the conditions given in Eqs, (11) and (12) 
and the constant radioactivity specified, i.e., 

[U]- [UJ, 0<f<r 1 

= [UJ, n<r< r 2 

= [UJ, r 2 < r < R 


The result is 


T(r,t)-T 3 =^V 


(h z + a ») sin <x n r 
Ra% + h(Rh + 1) 


[Ii(n) + I 2 (n)] exp ( — Ka%t) 


(15) 


where 

Ii(n) = Zii(n) + I„(») + h 3 (n) (16) 

lu(tt) — “T (sin a n r 1 — a n fi cos 4- |~ — ^ 4 — -1 (17) 

«« L r i a Tt J 


I 12 (n) 


[T. m (fl) + T' m R] 


(sin a ?i R — sin — ot n R COS a n R + a n r a COS «„fi) 


T 

v [2a n (R sin a n R — r t sin a n r J — (<4R 2 — 2 ) cos a n R 

<*n 

+ («» r i _ 2 ) COS «« r a)] ( 18 ) 
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(19) 


Z 13 (n) = f (sin a n R — a^R cos a n R) 

a n 


1 4 

J 2 (n) = — ~ Cj(fi,t) [S a i(/,fi) + S Z2 (/ ? rc) + S 33 (/,n)] 

^ i = i 


with 


( 20 ) 




1 — exp [ — (A; — Kc&)t ] 

(A / — KCtn)ctn 


( 21 ) 


S 21 (;,n) = ^!j(sin « H fi - cos a n fi) (22) 

S 22 (/,n) — 5 2 ^(sin a n r 2 — si — a n r 2 cos a w r 2 4- a n r x cos <x^r^) (23) 

S 2B (/,«) = Sa/sin a»R — sin «*r 2 — o^R cos + a n r 2 cos «„r 2 ) (24) 

The a n are consecutive positive roots of Eq. (10), The are given by 

Si i = pA {j H } exp (A,T), (/ = 1,2,3, 4) (25) 


where the A {j are abundances of species / in the ith shell and T = 4.6 X 10 9 years 
or an appropriate age for the object in question. 


The heat generation and decay constants used are taken from Clark (Ref. 8) 
and are listed in Table 1 in calories and in calories per gram units. 


Table 1. Radionuclide heat generation and decay constants 


; 

Nuclide 

Hj (cal seer 1 g -1 ) 

(seer 1 ) 

X 


6.658 X lO-^ 

1.682 X 10-i7 

2 

232 Th 

6.314 X 10-f 

0.1582 X 10~ 17 

3 

235U 

1.363 X 10-" 

3.082 X 10 17 

4 

23SU 

2.25 X 10- 8 

0.488 X 10" 17 


V. Solution for Change in Boundary Insulation at Time V 

If the insulation jchanges at time f > 0 to a value such that the constant h 
changes from h to h x (corresponding to insulation thickness changes d to d ly all 
other parameters remaining invariant), the solution for the subsequent tempera- 
ture is T(r>t)(t > f): 


9 

T(r,t) — r, = — ^ i/^n) si 


sin a n r 



X 


Ii(m) exp [— K(aW + ot«01 
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A 4 

+ ^ 2 ) *0 Q(1,m) 

+ V C'(/,n;<) g'(/,n) | 

1=1 > 




where 


*0») = 


ft 2 + ttm 

Ro& + /i(Rh + 1) 


(atn roots of a cot aR + h = 0) 


1,2 -U __ 

A (n) = =r , (an roots of a COt aR + h x = 0) 

Ral + MR/li + 1) 

h(m) = In(m) + J 12 (m) 4- l n (m) 


( 26 ) 


(27) 


(28) 

(29) 


where I n , I 12 , Z 13 , are given by Eqs. (17), (18) and (19), changing m for n. Note 
however that in Eq. (17), T m (fi) refers to the melting temperature at r,; there is no 
summation over m implied. Also in Eq. (26) 


C(?,m,n ; #/) = 


exp [ — x(«lf + alt)] — exp [~(Ajt' + 
(Xj — Kali) ali 


(30) 


C'(j,n-t) 


exp ( — — exp (—\jt) 

(A / — Ka%) a% 


(31) 


U(m,n) 


sin (am — <*»)R _ sin (»,» + «»)R 

(<*m — a„) (a m + or„) 


(32) 


Q(j,m) = S 21 (;,m) + S 22 (;,m) 4- S 23 (j,m) (33) 

S 21 (f,m) = ^ (sin a m r x cos atmfi) (34) 

S 22 (/,m) = s 2 >(sin atmf 2 — sin o^fi — a,„r 2 cos a m r 2 + o^fi cos a»,fi) (35) 

Sisd.m) — s 3 )(sin o^R — sin « m r 2 — o^R cos o^R 4- « m r 2 cos amf 2 ) (36) 

with the% given by Eq. (25). 
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The [Aj,] (i = 1,2,3 for core, mantle, and crust; j = 1,2, 3, 4 for nuclide species) are 
'U9Xl<Hj8[UJ 1.19X1(H/9[U,] 1.19 X 10* 4 ^[U 8 ] 


[An] — 


3.7 [U,] 

[Uil 

138.7 

137.7 


138.7 


Al2 

Aia 

a 14 


[U,] 

A Z i 

A 2 2 

Aj 3 

A Z i 


3.7 [U 2 ] 

[U g ] 

138.7 

137.7 


138.7 


[UJ 


A 32 

A 33 

A 34 


3.7 [U 3 ] 

[U 3 ] 

138.7 

137.7 


138.7 


[Ua] 


(37) 


and the relations between [U t ], [U 3 ], [U 3 ] as given by Eqs. (13) and (14) are 

Q'(i,n) = SJ, (/,«;*') + S »(/>;«') + S' 13 (f,n;t') (38) 


Sti(j,n;t') = sly(sin a„fi - a n r t cos a n r,) \ 

Sm,n;t f ) - s' y (sin — sin« n r t - a n r 2 cos a n r 2 + a n r x cos «„fi) \ (39) 

S 2 s(/,tX;#') = S3 j(sin a„R — sin 'X n r 2 — OnR COS a„R + a n r 2 COS a n r 2 ) \ 

s'ij = pAijHj exp [Ay(f - f')] (40) 

A "complete” discussion of a planetary thermal history thus involves a two-stage 
calculation with the present program: 1) the interval 0 < t < if, and 2) the interval 
f <t <T (present). 

Planetary density and [K]/[U] ratio are to be specified internally in the main 
program. 


VI. Input Data 

The following input data are required to make a calculation, in the order and 
format specified. 

(1) First card: Maximum number of time steps, maximum number of radial 

steps (1615). 

(2) Second card: Time increment A t (seconds); 2%, T 0 , T m (R), T' m (5D15.5). 

(3) Third card: r u r s , R, k, K (5D15.5). 

(4) Fourth card: Maximum value of ordinate for temperature in degrees Kelvin 

in plots (example, 3500°K) (5D15.5). 
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(5) Fifth card: d, K', t', d' (5D15.5). 

(6) Sixth card: [U^, /(5D15.5). 

The output consists of a printout of the input data (Fig. 3) and the temperature 
at specific times given as a function of radius (Fig. 4), and plots (Fig. 5), which 
each give initial temperature, as well as temperature at time t. 

The maximum number of radial steps in a calculation as well as the maximum 
number of iterations MAXIT in the main program (solutions to transcendental 
equations), are specified only in the main program. All input data are in calories 
and cgs units. Temperatures are printed in degrees Kelvin and times in seconds 
and years. Radii are given in centimeters. 

While we have dealt with homogeneous spheres, the solution given by Lowan 
(Ref. 2) is general enough to allow K, the thermal conductivity, to vary with 
radius. Whether the problem dealt with can be solved analytically depends upon 
whether expressions given in the original paper can be integrated. There is also 
no necessary restriction on forms for distribution of radioactivity; uniform istri- 
buttons have been used here lacking any reason to suppose otherwise in the Moon 
but exponential distributions could be handled as well. Problems involving thermal 
properties varying with temperature must be treated numerically. 

A listing of the main program and subroutines is given in the Appendix. 
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N 

3 


n 

20 


^"isJoaoo^ SURrACE I!"!!! "*!“!!* ''INITIAL* TEMPERATURE or ACCRETING MATERIAL MELT I NG POINT AT ZERO PRESSURE 

.1500000*017 .2500000*003 . <1500000*003 . 13*0000*P0<* 

MELTING POINT GRADIENT TMRl 'CORE* RADIUS 'MANTLE* RADIUS PLANETARY RADIUS THERMAL DIFFUSIVITY 

.Naoooou-oos » 1 R<*000fJ*0O'r . I <t b fl0 0 0*0 0 9 . 1 7 1 -000*00 * .I73SOOD-OOT ™e.M^D|FFUS I VITT 

THERMAL CONDUCTIVITY M A A I MUM ABSCISSA VALUF IN PLOTS 
.1100000*001 . 3r000Q0*00*( 

thickness of insulating later thermal CONDUCTIVITY OF INSULATING layer 
. 1 GO0U0O*0O* . I 000000*00- 

•T-PRJHE' new thickness OF INSULATING LATER 
•30QGOOO*017 ♦ J ODnCJQn + ou** 

ui r 

• 30QUOUO-00 7 • \ OOUODO^QO 1 


SOURCE terms 

. 1 7 7-43-0 I *4 

•ooooo 

•16174-013 


.29500-014 

•ooooo 
• 268?n-o 1 3 


•82040-01 4 

•aoooo 

• 74783-0 I 3 


• 4508S-Q J 4 
.00000 
*41097*013 


DECAY CONSTANTS 
K40: .1682000-016 

TN232 » .1882000-017 

U238: • 3082000-0 1 6 

U238I .4880000-01 7 

ROOTS OF EQUATION: TANULPH*R3) ♦ (ILPH/Hl * 0. 

1 . 1 07672323506-007 

2 *278012319980-007 

3 *456100055506-007 

4 .6356745807 11-007 

5 *815764264820-007 


Fig. 3. Listing of “input data" and first five roots of transcendental equation tan (a- R) + a /h = 0 



TEMPERATURE, °K 


TIME (SEC) 

TIME (YEARS) 

M0. OF 

.7500000+017 

.2378234+010 

TERMS IN SUMMATION 

RADIUS 

Tl+VS 


.8690000+007 

,2339513+04 

62 

.1738000+008 

.2341481+04 

47 

.2607000+008 

.2343264+04 

41 

.3476000+008 

.2344303+04 

38 

.4345000+008 

.2341946+04 

36 

.5214000+008 

.2333929+04 

28 

.6083000+008 

.2320373+04 

26 

.6952000+008 

.2292757+04 

26 

.7821000+008 

.2251677+04 

25 

.8690000+008 

.2193342+04 

25 

.9559000+008 

.2111307+04 

24 

.1042800+009 

.2001994+04 

24 

.1129700+009 

.1866203+04 

23 

.1216600+009 

.1700109+04 

25 

.1303500+009 

.1503824+04 

24 

.1390400+009 

.1279133+04 

27 

. 1477300+009 

.1026541+04 

27 

.1564200+009 

.7675109+03 

35 

.1651100+009 

.5251613+03 

41 

.1738000+009 

.2893931+03 

75 


Fig. 4. Listing of radius and temperature for accompanying plots for 
t = 7.5 X 10™ sec and [U] — 30 ppb 
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Appendix 
Program SPHERE 

«RUN»/R DEC 3 .J5N25Z. SPHERE .06.400/0. *183/516 . COREL 

• SC4020 BLOG/ 183 .BOX/516 • CAMERA/9 IN .FRAMES/9U 
■FOR. IS MAIM. MAIM 
C 

C PARAMETER* MAXR — MAXIMUM NUMBER OF RADIAL STEPS. 

C MAXIT — MAXIMUM NUMBER OF ITERATIONS 

C 

PARAMETER MAXR*40 

PARAMETER MAXIT-500 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION A(3>4) »H( A) 

COMMON /BLK1/ ALPHA (MAXIT) »ALPHAM< MAXIT ) »S ( 3 .4) »XLAMDA(4 ) . 

1 PSI (MAXIT) .PSI BAR (MAXIT) .XU (MAXIT). 

2 S2KMAXIT) .S22IMAX1T I .S23(MAXIT) . 

3 SBAR1 (MAX I T ) .SBAR21 MAXIT ) .SBAR3 (MAXIT I . 

4 TPR.T1ME.XKAPPA.XA0XX.R3.SPR (3.4) 

INTEGER JX(3).JY(3) »NP(3) .INTERPI 3) 

REAL XYIMAXR.4) >ROW2(2) 

REAL YMAX 

INTEGER SYMBOL ( 3) »T I TLE2 ( 14) .XNAME 1 14 ) . YNAMEilO I ,ROWl ( 2.2) 
DIMENSION I CONRG ( 4 ) .SUMO 14). SUMl 1 4 ) .KOUNT I 4 ) » IFLAGi 4 ) 

DATA XLAMDA/1. 682D-17.. 15820-17. 3. 0820-17.. 488D-17/ 

DATA PI/3. 14159265358979300/ 


DATA 

JX/1 *1 » 1/ JY/2* 

3*4/lNTERP/X*l*X/ 



DATA 

SYM6UL/6H1 

»6h2 *6H3 

/ 


DATA 

TITLE2/6H 

»6H »6H 

•6H1 . T1 

«6H »6H » 

1 

6H2. 1N1.6HTIAL T .6HEMPERA.6HTURE 

D»6HI STRIB»6HUT ION » 

2 

6H 

* 6H / 



DATA 

XNAME/6H 

* 6H *6H 

»6H 

t6H *6H » 

X 

6HRA01AL 

>6H DISTA»6HMCE 

*6H 

»6H *6H • 

2 

6H 

*6H / 



DATA 

YNAME/6H 

» 6H >6H 

$ 6H 

*6HTEMPERt6HATURE * 

1 

6H 

• 6H *6H 

>6H 

/ 

DATA 

IDIM/MAXR/ 





TAU* 1 . 435D1 7 
RH0=3.34D0 
BETA-2. D3 
H( 1 ) *6.6580-9 
H(2) *6.3410-9 
H13I-1. 3630-7 
H(4)»2. 2500-8 
ROW 1 < 1 » 1 1 *6HT I ME ( 

ROWl (2.1 )“6HSEC > 

R0W1 ( 1 .2 )*6HTIME ( 

ROW 1 (2.2) “6HYEARS ) 

EPS-l.D-3 

NCOUNT-9 

C 

C INPUT* N — NUMBER OF TIME STEPS 
C M — NUMBER OF RADIAL STEPS 

C 

10 READ! 5.500.END“400 ) N.M 
WRITE (6 .600 ) N.M 
IF(M.LE.MAXR) 00 TO 15 
MTEMP-MAXR 
WR1TE(6.607) MTEMP 
STOP 
C 

C INPUT* CELT — TIME INCREMENT 
C VS — SURFACE TEMPERATURE 

C TO — 'INITIAL* TEMPERATURE OF ACCRETING MATERIAL 

C TMO — MELTING POINT AT ZERO PRESSURE 

C TMPR -- MELTING POINT GRADIENT 

C R1 ~ ‘CORE* RADIUS 

C R2 -- ’MANTLE* RADIUS 

C R3 — PLANETARY RADIUS 

C X KAPPA — THERMAL DIFFUSIVITY 

C XK — THERMAL CONDUCTIVITY 

C ABMAX — MAXIMUM ORDINATE VALUE IN PLOTS (MINIMUM VALUE IS ASSUMEO 

C TO BE 0) 
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15 READ<5«501) DECT, VS»TO»TMu»TMPK ,Rl ,R2,R3 » X KAPPA ,XK 
REAl» (5* 501) ABMAX 

TMR1»TM0+TMPR*1R3-R1) 

WRITE (6. 601) DEL T » V S » TO » T MCI » TMPR * TMR 1 » R 1 » R 2 » R3 » XKAP PA » XX * ABMAX 
YMAX*ABMAX 

INPUT* 0 -- THICKNESS OF INSULATING LAYER 

XKPR — THERMAL CONDUCTIVITY OF insulating LAYER 
TPR — • T -PRIME* 

D1 — NEW THICKNESS OF INSULATING LAYER 

READ! 5*501 ) D»XKPR,TPR,Dl 
WRITE! 6 • 60S ) DtXKPR 
WRITE 1 6 , 610 1 TPR, 01 

INPUT* U1 
F 

READ ( 5,501 ) Ul»F 
WRITEI6,613) Ul,F 
U2*( l.DO-F)*Ul 

U3*U1*1 R3**3-R2**3+F* ( R2**3-R1**3 ) I / 1 R3**3-R2**3 ) 

FAC1*1.19D-A*BETA 

FAC2=137.7Du/138.7O0 

FAC 3*1. DO/ 136*700 

A(l»l>*FACl#Ul 

AI2»1)*FAC1*U2 

A( 3 * 1 ) *FAC1*U3 

At 1 ,2 ) ■3,7*U1 

A(2,2)=3«7*U2 

At 3,2 ) *3»7*U3 

A(1,3)-FAC3#U1 

A(2,3)»FAC3*U2 

A13»3)*FAC3*U3 

Atl*A)*FAC2*Ul 

A ( 2 , A ) *FAC2*U2 

A ( 3 , A ) =FAC2*U3 

DIF-TAU-TPR 

00 16 J*=l , A 

EX1*DEXP 1XLAMDA1 J I *TAU > 

EX2 * DEXP ( XL AMD A ( J t *01 F ) 

DO 16 I “1,3 

SU»J)«RH0*A|1»J)*H(J)*EX1 
SPR1 I,J)*RHO*Al I,J)*H1J)*EX2 

16 CONTINUE 

WRITE 16.602 ) 1 IS! I , J) »J“1*A) ,1*1,3) 

WRITE 16. 603) XLAMOA 
DO 17 1*1,3 
NP(I)*M 

17 CONTINUE 
XKOXK*XKAPPA/XK 
R1R1*R1*R1 
R3R3«R3*R3 
sigma»xkpr/ixk*di 

H0=SIGMA-1./R3 

HH*H0*H0 

TIMEsO.DO 

STATEMENT 999 DEFINES NEW TIME ORIGIN .T-DELT 
999 TIME*1.38758AD17 
0ELTR=R3/M 

OBTAIN ROOTS OF* TANIX*R3)+IX/HU)»0. 

00 19 1*1, MAX IT 

ALPHA l I ) *ALPH< 1 ,H0, R3 , JFLAG ) 

IF1JFLAG.EQ.0) GO TO 18 
WRITE16,609) 

GO TO 10 

18 AK*ALPHA 1 1 1 
AKAN*AK*AK 
AKR1=AK*R1 
«KR2*AK«R2 
AKR3*AK*R3 
AAR1R1*AKAK*R1R1 
SAKR1=DSIN(AKR1) 
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CAKR1-DC0SIAKR1) 

SAKR2=D$IN(AkR2) 

CALR2»DC0S( AKR2 ) 

SAKR3=DS 1 N( AKR3 I 
CAXR3=DC0S( AKR3 I 

PS1 ( I)*(HH+AKAK|/(R3*AICAM-H0*<R3*H0+1. I ) 

XI 11* ( TO/AKAK) *( SAKR1~AKR1*CAAR1 ) 

1 + ( ( TMR1-TO ) / ( AARlRl*AXAM )*<3.*t AAR1R1-2. )*SAKftl- 

2 AKRl*(AARlRl-6. )«CAXR1I 
XJ12={ CTM0+TMPR*R3)/AKAXI*ISAK.K3-SAK.K1-AK.R3*CAXR3+AK.R1*CAKR1) 

1 -(TMPR/(AX*AXAX| )*(2.*AMMR3*SAKR3-R1*SAXR1) 

2 -tAXAX*R3R3-2. l*CAXR3+(AARlftl-2. )*CAXR1) 
XI 13= ( -VS/AKAX 1 *( S AXR3-AXR3*CAXR3 1 

XI1(I)*XI11+XI 12+X 113 
S21< I »=SAXR1-AXR1*CAXR1 

S22I I ) = S A K.R2-SA K.R 1 “■ AXR 2*CA XR2 + AKR 1 *C AXR 1 
S23 ( I) =SAXR3-3AXK2-AXR3*CAKR3+AKR2*CAXR2 
19 CONTINUE 
190 SIGMA1*XKPR/(XK*D1| 

H1*SIGMA1-1./R3 

H1H1“H1*H1 

C 

C OSTAIN ROOTS OF' TAN < X*R3 >+<X/Hl 1 *0. 

C 

JO 230 I=1»MAX1T 
AUPHAMI n«ALPHII ,H1 .R3 ,JFLAGI 
1F( JFLAG.EO.OI 00 TO 210 
WRITE(6»609 ) 

60 TO 10 

210 AX»ALPHAM( I ) 

AXAX=AX*AX 

AKR1»AK*R1 

AXR2-AK*R2 

AKR3=AX*R3 

SAKR1«DSIN(AKR1) 

CAKR1*DC0S( AKR1 > 

SAKR 2*05 1 N ( AKR2 ) 

CAKR2*DC0S I AKR2 ) 

SAKR3=DS I N ( AKR3 > 

CAXR3-DCOSIAXR3) 

PS1BARC I )*(H1H1+AKAK)/(R3*AKAX+H1*(R3*H1+1 . > I 
SdARlII )=SAXR1-AXR1*CAXR1 

SOAR 2 < 1 1 »0AXR2-SAXR1-AXR2*CAXR2+AXR1*CAXR1 
SOAR 3(11 =SAXR3-SAKR2-AK,R3*CAKR3+AKR2*CAXR2 
230 CONTINUE 

00 280 1*1, N 'LOOP OVER TIME 

T1ME-T1ME+DELT 

TIM£YR=TIME/3.1536D07 

R0W2( 11= TIME 

R0N2C 2) -TIME YR 

R*»0*DU 

WRi T£ I 6*611 ) TIME.TIHEYR 

00 270 J* 1 »M • LOOP OVER RADIAL DISTANCE 

r»r+oeltr 

c 

C INITIALIZE convergence flags and sums 

c 

1C0NRG( 1 )“0 
KOUNTUjttO 

1 FLAG ( 1 ) *0 
SLIM0(1)»0.D0 
SUHI ( 1)*»0.00 

DO 250 K=1*MAXIT 
\K=ALPHAHIK) 

AKA<*AK*AK 

AXR«AX*R 

SAJCR°05IN(AKR> 

Xl>XlEX(AAAK f AK»K) 

SUM1 ( 1)*SUM1(1 )+PSIdAR(K)*SAICR*Xl 

IF(DABS( (SUMl< 1 )-$UMO< ln/SUMl< 1 n *GT*EPS) GO TO 240 

IF{ IFLAGU) .NE.KM1J KOUNTll)*l 

IF( I FLAG < 1 ) *£u* kN I ) KOUNTIl > * KOUNTll ) +1 

IFLAGtlMK 

240 SUM0m»SUMlU) 

ICONRGU)«K 

IF( X.OUNT ( 1 ) *GT gNCOUNT ) GO TO 260 
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250 CONTINUE 
260 T1“2.*SUM1I 1>/R 

XY! J*1)“R 
XY ! J»2 ) “Tl+VS 
IF(R.GT.Rl) GO TO 262 
XY(J»3)“TQ+( TMftl-TO J*R*R/R1R1 
GO TO 266 

262 XY( J»3I“TM0+TMPR*(R3-R) 

266 WR1TE16.612) R »X Y ( J .2 ) • 1C0NRG (1 ) 

270 CONTINUE , , „ 

CALC KCPL I XY » IDIM* 0 »UX»UY *NP » 1NTERP .SYMBOL * TITLE2 »XNAME» YNAME » 

X ROW I » ROW 2 » 2 t 

CALL GR10Y!-.1»YMAX»0,0.0.6H 
CALL KCPL1!2,2»1} 

260 :ONTINUE 
GO TO 10 
600 STOP 

500 FORMAT (161 5) 

501 FORMAT (5015. 5) 

502 FORMAT ( 6015*5) 

600 FORMAT! *0 N M'/1H .215) 

601 FORMAT ( • OTIME INCREMENT SURFACE TEMPERATURE • 'INITIAL" TEMPER 
1ATURE OF ACCRETING MATERIAL MELTING POINT AT 2ER0 PRESSURE'/IH i 
2D14.7,8X,D14.7,32X»D14.7.19X.014.7/'0MELTING POINT GRADIENT 

3 TMR1 "CORE" RADIUS "MANTLE" RADIUS PLANETARY RADIO 

AS THERMAL OlFFUSl TIVJTY'/IH .8X.D16.7.2X*D16.7»3X.D16.7»6X»D16.7 
5,5X.Di6.7»lLX.D16.7/'0THERMAL CONDUCTIVITY MAXIMUM ABSISSA VALUE 
6 IN PLOTS'/IH ,6X,D14.7.19X.D14.7> 

602 FORMAT ( ' OSOURCE TERMS'/ 1H ,4015. 5/1H *6015. 5/1H .6D15.5) 

603 FORMAT! '0 DECAY CONSTANTS'/' K40 " »4X»D14. 7/ • TH232"*2X» 

1 016.7/ • U235" .3X.D16.7/' U238 " .3X.014. 7 > 

606 FORMAT! '0 TIME'/IH *D16.7/'0 RADIUS Tl+VS 

1 TR+VS V“T I +TR+V5' ) 

605 FORMAT (1H »4D14.7 ,41 10 ) „ , 

606 FORMAT! '0 PARTIAL! TI) PARTIAL! TR) FLUX 1 FLUX 2 

1 TOTAL FLUX' /1H ,5014.71 

607 FORMAT (' 0NUM6EK OF RADIAL STEPS EXCEEDS MAXIMUM DIMENSION'/ 

1 • MAXIMUM DIMENSION" ,141 

60S FORMAT ( • OTHICKNESS OF INSULATING LAYER THERMAL CONDUCTIVITY OF I 
INSULATING LAYER'/IH »15X.D14.7,29X,D14.7J 

609 FORMAT OOFAULTY ROOT OBTAINED FROM SOLUTION TO' TAN ( X*R3 1 + (X/R3 I *0 
1. COMPUTATIONS TERMINATED. • I 

610 FORMAT! *0 "T-PRIME" NEW THICKNESS OF INSULATING LAYER'/ 

1 1H »D14.7»23X,D14.7I 

611 FORMAT! '0 TIME (SEC) TIME (YEARSI'/IH , D14. 7.5X ,014.7/ 

1 1H0 ■ SX, ' RADIUS' »9X » 1 Tl+VS ' I 

612 FORMAT 1 1H ,D14.7,E14.7,28X,110) 

613 FORMAT ( * 0 U1 F'/1H ,2D14.7) 

END 

•FOR, IS X1EX.X1EX 

DOUBLE PRECISION FUNCTION XlEX(AKAk,AK,K) 

PARAMETER MAXlT«500 

IMPLICIT DOUBLE PRECISION (A-H.O-2) 

COMMON /BLK.1/ ALPHA (MAX I T I , ALPHAM! MAXI T I ,S 1 3 ,4) ,XLAMDA !4 ) » 

1 PS1 1 MAXI T ) , PS I BAR! MAX I T ) ,XI 1 (MAXI T ) » 

2 621 1 MAXI T > .S22 (MAXI T > ,S23 < MAX l T ) . 

3 SBARUMAXIT! .SBAR2IMAXIT I »SBAR3 (MAX IT > » 

4 TPR .TIME .XKAPPA ,XKOXK,R3 ,SPR (3,4) 

DATA EPS/1. O-S/NCQuNT/9/ 

KOUNT«0 

IFLAG“Q 

SUH0“v.D0 

SUM1»U.D0 

DO 20 I«1,MAXIT 

AI “ALPHA (1 1 

AIAI“AI#A1 

0IF-A1-AK 

SUM“AI+AK 

SDR3“DSIN1DIF*R3> 

SSR3“DSlN(SUM*R3t 

POWcRl“-XKAPPA*AIAI*TPR-XXAPPA*AXAA*TIME 
P0WER2“— XKAPPA* ( A1 AI *TPR+AKAK*T IME ) 

EX«DEXP!P0WER1> 

EX2“0EXP ( P0WER2 ) 

XI1EX«X11(1I*EX 

XI2EX“0.D0 
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00 8 J*1*A 

P0WER3*- * XLAMDAt J ) *TPR+XKAPPA*AKAK*T IME 1 

CBA*U < EX2-0EXP ( PQWER3 ) ) / < < XLAMOA * J ) -XKAPPA*Al AI )*AI A l> 

Q-$*1»J)*S21< I)+S*2*J)*S22< 1 )+S<3# J >*S23 * i ) 

X12EX«XI2EX+C*AR*0 
a :qntinue 

XI2EX«XI2EX*XKOXK 

SUM1*SUM1+PSK I)**XI1EX+XI2£X)** ISDR3/DIF > -*SSR3/SUM > ) 

IF l DABS* t SUMl-SuMO ) /SUM1 ) *GT .EPS ) GO TO 10 
IMIM-l 

lFUFLAG.NEilMl) KOUNT-1 
I F < 1 FLAG«EQ« IM1 j KOUNT*KQUNT+ 1 

1 FLAG* I 

I F ( KQUNT *GT •NCGUNT I GO TO AO 
10 SUM0*SUM1 
20 CONTINUE 

WRITE * 6 » 30 ) 

30 FORMAT ( * GFAiLURE TO CONVERGE IN X ibARl ' ) 

AO P0w£R2“-XKAPPA*AKAK*TIME 
EX2»DEXP<PQ*ER2) 

DO 50 J*1 f A 

POWER3— XlAMDA(J)*TIME 

CbARPR* ( EX2-0EXP ( POWER3 ) ) / * * XLAMOA * J ) -XKAPPA*AKAK J *AKAK ) 
QPR*$PR * 1 • J )*SBAR1 ( X )+SPR t 2 p J )*5BAR2 t X I +SPR * 3fU )*SBAR3 (K) 
SUM1»SUMH-XKQXK#CBARPR*GPR 
50 CONTINUE 
X1EX-SUM1 
RETURN 
ENO 

1 FOR 1 1 S ALPH* ALPH 

DOUBLE PRECISION FUNCTION ALPH U ,H .R3 , JFLAG > 

IMPLICIT DOUBLE PRECISION (A-HtQ-Z) 

FUNCTION ROUTINE WHICH CALCULATES SUCCESSIVE POSITIVE ROOTS OF THE 
FUNCTION* TAN< ALPHA*R3 ) + * ALPHA/H) * 0, 

DATA PI / 3 * 1 A159 265 3 589 7 9 32 ADO / 

DATA N COUNT/A/ 

kount«i 

lChECK*0 

IF(X.NE.I) GO TO 10 
WRITE (6* 5) 

5 FORMAT* ‘OROOTS OF EQUATION* TAN*ALPH*K3) + (AlPH/H) * 0* » ) 
JFLAG=0 
STEP=PI/R3 
STEP02*STEP/2# 

EPS«H/1000. 

X=STEPOZ+EPS 

oneoh=u/h 

GO TO 20 

10 X*K*STEP-STEP02+EPS 
20 DO 30 I«l»200 
IM1-I-1 
THET A c X*R 3 
TEMP=l./DCOS* THETA) 

F-DTAN*THETA)+X/H 
FPR“R3*T £MP*TEMP+QNEGH 
XNEW-X-F/FPR 

IF* DABS* *XNEW“X) /XNEW l *GT* 1 *0-12 ) GO TO 25 
IF* ICHEOUNE.IMI) KOUNT=l 
IF * I CHECK* EG* I MI ) KOUNT-KOUNT+1 
I CHECK- I 

I F ( KOUNT *GT *NC0UNT ) GO TO 50 
25 X-XNEW 
30 CONTINUE 

WRITE* 6 >40 ) XtXNEW.F 

AO FORMAT* 'ONEaTON METHOD FAILED TO SATISFY TERMINATING CRITERIA 
I X XNEW F*/1H »52X»3D16»8) 

STOP 

50 ALPH-XNEW 

WRI TE *6» 55) K >XNEW 
55 FORMAT* 1H »Il0iD22tl2) 

%MID*K*P I/R3 

XLBND*XMID-STEP02 

UBOUND»XMID+STEP02 

IF(XLBND«LT*X*AND*XoLToUBGUND> RETURN 


JPL TECHNICAL MEMORANDUM 33-718 


ORIGINAL PAGE, IS 
OF POOR QUALITY 



WRITE<6*60) K»XL0ND tXtUBOUND 

60 FORMAT! *0S0LUT ION NOT PROPERLY BOUNDED K LOWER BOUND 

1 X UPPER BOUND •/ lh »3QX » 1 5 » 3D16 « 6 ) 

JFLAG*1 

RETURN 

END 

♦ MAP » L 

Lib Llb*PLOTS 

• XQT 


3 20 


1.5016 

230* 

450 • 


1360* 

4*80-6 

1. 48600 S 

1*7 14008 


1 .736008 

1*640-2 

1*100-2 

350u . 






100000. 

•00001 


* 3D+17 

1000* 


.30000-07 

1 • 





1.5016 

130* 

250* 


1373 

5*00—6 

3 20 






1.5016 

230* 

450* 


1360. 

4.80-6 

1.486008 

1*714008 


1 *738006 

1*640-2 

1*100-2 

3500. 






100000. 

•00001 


*30+17 

1000* 


.4000D-07 

1. 
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